The DBD configuration is depicted in Fig. 1. A high voltage electrode is located at y = 0 mm and a grounded electrode is located at y = 10 mm, with a dielectric layer covering on its surface (the thickness of the dielectric layer is assumed to be zero). Thus, the discharge gap has a size of d_G = 10 mm. The discharge varies only along the y direction to establish a one-dimensional (1D) model (the discharge is uniform in the x direction as we ignore the edge effect and others). The gap is filled with a premixed stoichiometric hydrogen–oxygen mixture at an initial pressure of 1.0 atm (1 atm = 1.01325 × 10⁵ Pa), and at a temperature of 500 K.

Fig. 1.

	Figure Option ViewDownloadNew Window
	Fig. 1. (color online) Schematic diagram of discharge simulation configuration.

The governing equations for discharge simulation consist of the drift-diffusion equations for positively charged, negatively charged, and neutral particles, and Poisson’s equation for the electric field as follows: 1 2 3 4 where n_k is the number density of particle k; μ_k and D_k denote the mobility and diffusion coefficients respectively; S_k is the reaction source term; φ, e, ε₀, and ε_d are the electric potential, elementary charge, permittivity of free space, and relative permittivity, respectively; ε_d = 3.0; n_k and Z_k are the number density and charge of particle k, respectively. Then, the strength of the electric field is calculated by the following equation: 5 A total of 47 discharge elemental reactions^[25] are used where the positively charged particles O⁺, , , H⁺, negatively charged particles O⁻, , , OH⁻, and neutral particles O, O₂, O₃, H₂, H, OH, H₂O are taken into account. Boundary conditions are expressed as n_k = 0 (for negative charged particles) and ∂n_k/∂y = 0 (for positive charged particles) at the high voltage electrode, ϕ = 0 V at the grounded electrode, and ∂n_k/∂y = 0 at the surface of the dielectric layer. The particle velocity at the electrode surface is limited in practice, yet the boundary conditions used here lead to an infinite velocity value of particle, so the distributions of species concentration are given by fitting functions, which are determined based on the overall distributions excluding the several values near the electrode boundary.

The DBD actuator is driven by a sinusoidal AC power supply and the frequencies of 8 kHz, 10 kHz, and 12 kHz are selected for comparison, while the voltage is fixed at 14 kV. The selected actuating frequencies are in the range of the commonly used frequencies in our AC-DBD experiments. A quasi-neutral plasma cloud is used as the initial condition; the number densities of and electrons are both 10¹⁵ m⁻³ initially in the entire domain, while other species are assumed to be zero.

Poisson’s equation is discretized by the central differential method and the successive over relaxation (SOR) approach is utilized. The Crank Nicolson-finite element method (CN-FEM) is used to solve the drift-diffusion equations. The convective and diffusion terms are evaluated by the upstream scheme and central differential scheme, respectively. The spatial and temporal computational grid sizes are 5 μm and 0.02 ns, respectively. Details concerning the equations, reaction mechanisms, boundary and initial conditions, together with the numerical solving schemes can be found in Ref. [25].

Based on the obtained discharge results, an arrangement of electrode, which is compatible to the discharge simulation in the first stage, is put forward for the detonation tube as shown in Fig. 2. The entire computational domain of reaction flow is a cylinder tube with an axial length of 800 mm. The inner diameter is 23 mm with a discharge gap size of 10 mm. The closed side wall, also called the thrust wall of a detonation combustor, is on the left end of the tube and the outlet is at the right end. A rod shaped high voltage electrode is connected to the center of the closed side wall with an axial length of 10 mm and an outer diameter of 3 mm. The cylindrical wall of the detonation tube is set to be the grounded electrode directly with a dielectric layer covering its internal surface, and the axial length of the grounded electrode is also 10 mm. To restrict the discharge region between the two electrodes (this domain is called “the discharge space”), insulating material is adopted for the other part of the wall. For simplicity, the discharge is assumed to be uniform along the axial and circumferential directions with neglecting the influence of curvature on the spatial electric field, which is generally acceptable as reported in Refs. [26] and [27], so the former discharge results can be used in the later simulation. For easy comparison, the computational configuration stays unchanged when the DBD actuator is turned off, and the traditional heat ignition is performed within the discharge space for all cases.

Fig. 2.

	Figure Option ViewDownloadNew Window
	Fig. 2. (color online) Computational domain and electrode arrangement.

The loosely-coupled approach is used to deal with the simulation of plasma assisted detonation initiation here.^[25,28] According to the operation sequence, that is, the first step is to start the discharge and the second step is to perform the heat ignition, the spatial distributions of key radicals and molecules at the moment when all the radicals reach their relatively high concentrations, which is obtained in the first stage (i.e., discharge simulation), together with the widely used traditional heat ignition,^[29] are chosen as the initial conditions of the second stage (i.e., combustion simulation) in the discharge space.

Structural meshes are adopted in the whole computational domain. One mesh scheme that consists of 839 800 cells is chosen, with the minimum grid step being 0.1 mm. An axisymmetric boundary is specified for the axis of the detonation tube, so the total mesh cell number is half of the above-mentioned original mesh cell number. When the outflow is subsonic, the local static pressure is set to be 1 atm at the outlet and all other variables are extrapolated from the upstream zone; when the outflow is supersonic, a supersonic extrapolation condition is assumed here. The species is assumed to be air at the boundary. No slip and adiabatic conditions are imposed on the entire wall of the tube. The tube is filled with stoichiometrically mixed quiescent hydrogen and oxygen with an initial temperature of 500 K and a pressure of 1.0 atm. The heat ignition zone, which is adjacent to the closed side wall, has a temperature of 3000 K and an axial length of 10 mm.

The unsteady Reynolds averaged Navier–Stokes equations together with the realizable k–ε turbulence model, chemical source term, and multi-species conservation equations are solved as reaction flow governing equations. The finite rate reaction model is used for combustion simulation, and the 6-species 7-reaction scheme is selected as the O₂–H₂ chemical kinetic model. The dynamic processes of detonation wave formation and its propagation are both well captured by the simulation program, and details about related equations, physical models, solving methods, and validation can be found in Refs. [25] and [30].

In order to study the influences of the actuating frequency on the temporal and spatial distribution of discharge plasma, three frequencies including f = 8 kHz (Case 1), f = 10 kHz (Case 2), and f = 12 kHz (Case 3) are chosen. The temporal and spatial distributions of the number density of key active radicals O, H, and electrons in the 4th discharge cycle T for the three cases are shown in Fig. 3. It is found that the frequency does not change the pattern of the distribution of every particle basically, yet it changes the magnitude of number density. The density of every particle declines as the frequency rises from 8 kHz to 12 kHz. This should be attributed to the fact that a higher actuating frequency means a shorter duration of single discharge, leading to a decrease in the ionization degree. Furthermore, atom O is the most sensitive to the variation of frequency according to the relative change of the peak value of particle density. Specifically, the maximum density of Case 3 is reduced by about 43.3% compared with Case 1. As for the other two particles, it is found that only the electron density decreases in constant amplitude as the frequency rises continuously.

Fig. 3.

	Figure Option ViewDownloadNew Window
	Fig. 3. (color online) Temporal and spatial distributions of discharge products. (a) f = 8 kHz, H, (b) f = 10 kHz, H, (c) f = 12 kHz, H, (d) f = 8 kHz, O, (e) f = 10 kHz, O, (f) f = 12 kHz, O, (g) f = 8 kHz, electron, (h) f = 10 kHz, electron, (i) f = 12 kHz, electron.

The reason for the concentration of active species decreasing as theactuating frequency rises is related to the discharge power, which can bedetermined by the following equation: 6 where U_i and I_i are the voltage and the current at the i-th moment, n is the total number of moments of one cycle.

As shown in Table 1, the discharge power decreases as frequency rises in a nonlinear manner. When other discharge environmental and actuator control parameters are kept unchanged, the energy input into the plasma will increase with a larger discharge power. This will benefit the break of the species chemical locks, such as chemical bonds, then more active particles are formed.

Table 1.

Table 1.

Table 1. Relationship between actuating frequency and discharge power. . f/kHz 8 10 12 P/W 14.9 12.8 11.1

Table 1. Relationship between actuating frequency and discharge power. .

Among all the intermediate species of the discharge plasma, radicals O, H, and OH have notably higher number densities than the others in the simulation results (other intermediate species are not presented here). It has been reported^[1,31] that it is the three species that have relatively high concentrations in a hydrogen-oxygen discharge plasma, and have the strongest influence on ignition under the given conditions. Therefore, only species O, H, and OH are considered in the combustion simulation. In view of the fact that nearly all of the three key radicals reach their highest concentration at about t = 0.6T based on the temporal and spatial distributions of their density for every case, the concentration distributions of O, H, OH, H₂, and H₂O at 0.6T are therefore selected as the initial conditions of subsequent combustion simulations, while other species with low concentration are ignored except for O₂.

As is well known, the whole DDT process consists of five sub-processes, which are of the very slow combustion right after ignition, deflagration with a rising burn velocity, over-driven detonation, decay of over-driven detonation, and propagation of self-sustain stable detonation. Three flow statuses including the static, subsonic, and supersonic exist in the detonation tube during a DDT history. Thus, both the flow field at typical moments and the whole history of flow evolution are required for comparison and analysis.

3.2.1 Before the formation of stable detonation waves

Figure 4 shows the distributions of Mach numbers in the detonation tube at two typical moments t = 70 μs and t = 100 μs. These two moments are selected because they reflect the flow behaviors at the subsonic stage and the supersonic stage, respectively, before a stable detonation wave forms, and it is easy to distinguish the wave front from a Mach number distribution at the moments. The flow field structures of the three cases which have different actuating frequencies are almost the same at a certain moment. The flow is still in subsonic state at t = 70 μs yet it reaches supersonic at t = 100 μs. It is found that the wave front becomes close to the thrust wall for the given moments as the frequency rises from 8 kHz to 12 kHz by comparing the x positions of the leading wave in the panels, indicating that the DDT process develops slower as the frequency rises, however the change is very small.

Fig. 4.

	Figure Option ViewDownloadNew Window
	Fig. 4. (color online) Mach number distributions in detonation tube under different frequencies. (a) f = 8 kHz, t = 70 μs, (b) f = 8 kHz, t = 100 μs, (c) f = 10 kHz, t = 70 μs, (d) f = 10 kHz, t = 100 μs, (e) f = 12 kHz, t = 70 μs, (f) f = 12 kHz, t = 100 μs.

To realize the accumulated variation led by the plasma during the whole subsonic phase, the profiles of pressure, temperature, and species concentration along the tube axis at Ma = 1 are presented in Fig. 5. An x range covering the abrupt changing curves of flow field parameters is selected here. Like the description of the flow structures in Fig. 4, none of the overall profile characteristics of each parameter are changed by the frequency at Ma = 1. Only the positions of the leading shock waves of the three cases are different. Specifically, the leading shock wave in the f = 12 kHz case moves farthest away from the thrust wall, the next is the f = 10 kHz case, and the last is the f = 8 kHz case. The fine differences can be distinguished from Fig. 5(c) clearly: the distances in x-direction between every two adjacent steep curves, which represent the reaction boundary, are both about 0.4 mm.

Fig. 5.

	Figure Option ViewDownloadNew Window
	Fig. 5. (color online) Profiles of flow field parameters along the x axis at Ma = 1 under different frequencies: (a) pressure, (b) temperature, (c) mass fraction of H2 and H2O.

3.2.2 DDT time and distance, and steady propagation stage

Figure 6 shows plots of DDT time and distance versus actuating frequency, which are received from the profile of the flow field parameter along the x axis at the moment when the Chapman–Jouguet (CJ) detonation wave has just formed. It is clear that the time and distance required for the stable detonation wave forming increase as frequency rises. The change caused by the frequency rising from 10 kHz to 12 kHz is notably larger than that caused by the frequency rising from 8 kHz to 10 kHz. Specifically, the former is about 3.9 times the τ_DDT of the latter, and it is about 4.9 times the L_DDT of the latter, indicating that the DDT time and distance increase in a non-constant amplitude way, though the frequency rises in an equal amplitude.

Fig. 6.

	Figure Option ViewDownloadNew Window
	Fig. 6. (color online) Plot of DDT time and distance versus frequency.

To further realize the effect of actuating frequency on the propagation of a stable detonation wave, the profiles of pressure, axial velocity, and species concentration along the tube axis at t = 260 μs are shown in Fig. 7 for the three cases. At this moment, the abruptly changing curves of the flow parameters become steeper than those of Fig. 5, and the peak values of pressure profiles, especially the profiles of the mass fraction of H₂ and H₂O, are higher.

Fig. 7.

	Figure Option ViewDownloadNew Window
	Fig. 7. (color online) Profiles of flow field parameters along the x axis at t = 260 μs under different frequencies: (a) pressure, (b) axial velocity, (c) mass fraction of H2 and H2O.

In accordance with the detonation initiation orders (i.e., the DDT time and distance) of the three cases, the detonation wave of the lowest frequency (f = 8 kHz) case moves farthest at this moment, the next is the f = 10 kHz case, and the last is the f = 12 kHz case. In addition, the length relationship of the distance between two adjacent leading shock waves among the three cases also accords with the length relationship of the DDT distance, that is, the wave front of the f = 12 kHz case lags behind a lot, while the interval of the wave front between the f = 8 kHz case and the f = 10 kHz case is very narrow. Compared with the positional relationships among the wave fronts of these cases at the moment when a stable detonation wave has not yet formed (as shown in Fig. 5), the distance between the f = 12 kHz case and the other two cases is elongated notably, so the positional relationships of the wave front at the three frequencies are changed after the completion of the detonation initiation. In summary, the effect of the actuating frequency on the detonation forming process varies in the five sub-processes.

3.2.3 History of thrust wall pressure

Due to the close correlation between the thrust wall pressure and the performance of detonation engine, and the prominent unsteady characteristics of the pressure, it is necessary to study the dynamic characteristics of the pressure on the thrust wall.

Figure 8 shows the history of thrust wall pressure under different actuating frequencies. Here, the thrust wall refers to the inner closed side wall of the tube excluding the surface of the central electrode, and the magnitude of the pressure is an area weighted average value. Generally, the frequency seems to have no influence on the history. The pressure eventually tends to a steady value after a remarkable summit. According to the magnification of the highest pressure summit, the f = 8-kHz case reaches its maximum value first, the next is the f = 10-kHz case, and the last is the f = 12-kHz case. From the time span shown in the magnification graph, we can see that the flow field is still in the detonation evolution phase, and the axial differences among the three cases are nearly the same, which are in accordance with the above mentioned results about the influence of frequency on flow field evolution before a stable detonation wave is formed. As for the amplitudes of pressures, there is no difference among them. Thus, it is assumed that the actuating frequency does not affect the amplitude of thrust wall pressure.

Fig. 8.

	Figure Option ViewDownloadNew Window
	Fig. 8. (color online) Histories of pressure on thrust wall under different frequencies.

Based on the analysis above, it is concluded that a relatively low actuating frequency makes it faster to initiate the detonation, and the influence of frequency on the flow field evolution varies with time nonuniformly. However, the amplitudes of a certain flow field parameter are almost identical for different frequencies in the whole evolution process.